#####################
## The code in this file provides R code to
## reprpduce the analysis and graphics in #
## Ahlquist & Wibbels AJPS "Riding the Wave"
## Contact john.ahlquist@gmail.com with questions
## or problems.
## Last updated: 23 Nov 2011
#####################

library(foreign)
library(MASS)
library(WhatIf)
library(xtable)
library(car)
library(lmtest)
library(cshapes)
library(spdep)
library(RColorBrewer)
library(lme4)
library(Design)
library(sandwich)
library(R2HTML)
library(ROCR)
library(verification)



set.seed(98)

md<-read.csv("AhlquistWibbelsRidingTheWaveAJPSdata.csv")
md$prop.dem.svolik<-100*md$prop.dem.svolik

is.rs<-subset(md,country.existence==1 & year >1874 & WB.code != "ZAR") #ZAR has odd data for labor endow
noheg<-unique(md$WB.code)
is.rs.noheg<-subset(is.rs, WB.code != "GBR" & WB.code != "USA") #the "no hegemon" data frame 


## Figure 1
temp<-subset(is.rs,WB.code=="USA")
par(mfrow=c(2,2))
#pdf("filename.pdf")
plot(temp$year,temp$world.trade.gdp.interp,type="l", lwd=3, ylim=c(0,60),
     main="World trade & democracy (Polity IV)", xlab="year", ylab="percent")
lines(temp$year,temp$prop.dem.polity, lwd=3, col=grey(.5))

plot(temp$year,temp$world.trade.gdp.interp,type="l", lwd=3, ylim=c(0,60),
     main="World trade & democracy (Svolik)", xlab="year", ylab="percent")
lines(temp$year,temp$prop.dem.svolik, lwd=3, col=grey(.5))
dev.off()

## prep work for figures 2 & 3
region<-read.csv("./data to merge/WBcountryClass.csv",header=F) # uses external data file; be sure to adjust filepa

temp.65<-subset(is.rs,year==1965)
temp.00<-subset(is.rs,year==2000)
ins<-temp.65$WB.code[temp.65$WB.code%in%temp.00$WB.code]
forRLE<-data.frame(matrix(NA,nrow=length(ins),ncol=1))
rownames(forRLE)<-ins

temps<-subset(is.rs, WB.code%in%ins & year==1965, select="svolik.dem")
forRLE[,1]<-temps
temps<-subset(is.rs, WB.code%in%ins & year==2000, select="svolik.dem")
forRLE[,2]<-temps
temps<-subset(is.rs, WB.code%in%ins & year==1965, select="labor.endow")
forRLE[,3]<-temps
temps<-subset(is.rs, WB.code%in%ins & year==2000, select="labor.endow")
forRLE[,4]<-temps
forRLE[,5]<-region[match(rownames(forRLE),region[,1]),2]
temps<-subset(is.rs, WB.code%in%ins & year==1965, select="polity.democracy")
forRLE[,6]<-temps
temps<-subset(is.rs, WB.code%in%ins & year==2000, select="polity.democracy")
forRLE[,7]<-temps
colnames(forRLE)<-c("dem.65","dem.00","rle.65","rle.00","region","pol.65","pol.00")
forRLE<-na.omit(forRLE)

nclr<-length(unique(forRLE$region))
palette(brewer.pal(nclr,"Dark2"))

#Figure 2
#pdf("filename.pdf")
plot(log(forRLE$rle.65), log(forRLE$rle.00),col="white",xlab="1965", ylab="2000",main="RLE by region")
text(log(forRLE$rle.65), log(forRLE$rle.00),labels=rownames(forRLE), col=unclass(forRLE$region))
abline(h=mean(mean(log(forRLE$rle.00))),lty=2)
abline(v=mean(mean(log(forRLE$rle.65))),lty=2)
dev.off()

#Figure 3
#pdf("filename.pdf")
par(mfrow=c(2,2))
palette(c("red", "black"))
plot(log(forRLE$rle.65), log(forRLE$rle.00),col="white",xlab="1965", ylab="2000",main="RLE by 1965 democracy (Svolik)")
text(log(forRLE$rle.65), log(forRLE$rle.00),labels=rownames(forRLE), col=forRLE$dem.65+1)
abline(h=mean(mean(log(forRLE$rle.00))),lty=2)
abline(v=mean(mean(log(forRLE$rle.65))),lty=2)

plot(log(forRLE$rle.65), log(forRLE$rle.00),col="white",xlab="1965", ylab="2000",main="RLE by 2000 democracy (Svolik)")
text(log(forRLE$rle.65), log(forRLE$rle.00),labels=rownames(forRLE), col=forRLE$dem.00+1)
abline(h=mean(mean(log(forRLE$rle.00))),lty=2)
abline(v=mean(mean(log(forRLE$rle.65))),lty=2)

plot(log(forRLE$rle.65), log(forRLE$rle.00),col="white",xlab="1965", ylab="1965",main="RLE by 1965 democracy (polity)")
text(log(forRLE$rle.65), log(forRLE$rle.00),labels=rownames(forRLE), col=forRLE$pol.65+1)
abline(h=mean(mean(log(forRLE$rle.00))),lty=2)
abline(v=mean(mean(log(forRLE$rle.65))),lty=2)

plot(log(forRLE$rle.65), log(forRLE$rle.00),col="white",xlab="1965", ylab="2000",main="RLE by 2000 democracy (polity)")
text(log(forRLE$rle.65), log(forRLE$rle.00),labels=rownames(forRLE), col=forRLE$pol.00+1)
abline(h=mean(mean(log(forRLE$rle.00))),lty=2)
abline(v=mean(mean(log(forRLE$rle.65))),lty=2)
dev.off()

#Model 0.
## DV = Svolik
## omits RLE*Trade interaction for comparison purposes 
sv.ni.rle<-glm(svolik.dem ~ I(svolik.dem.m1) + log(labor.endow) + world.trade.gdp.interp +
              I(svolik.dem.m1*log(labor.endow))+ I(svolik.dem.m1*world.trade.gdp.interp) +
              prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) +
              nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) +
              prev.dem.death.sv + I(svolik.dem.m1*prev.dem.death.sv) +
              communist + as.factor(xr.reg) + I(nb.demtrans.sv.m1>0) + I(svolik.dem.m1*gdppc.growth),
              family=binomial(link="probit"), data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))

s2.sv.ni.rle<-vcovHC(sv.ni.rle)# "robust" SEs
AIC(sv.ni.rle,k=log(length(sv.ni.rle$residuals))) # 1406


#Model 1.
sv.i.rle<-glm(svolik.dem ~ I(svolik.dem.m1) + log(labor.endow)*world.trade.gdp.interp +
              I(svolik.dem.m1*log(labor.endow))*I(svolik.dem.m1*world.trade.gdp.interp) +
              prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) +
              nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) +
              prev.dem.death.sv + I(svolik.dem.m1*prev.dem.death.sv) +
              communist + as.factor(xr.reg) + I(nb.demtrans.sv.m1>0) + I(svolik.dem.m1*gdppc.growth),
              family=binomial(link="probit"), data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))
AIC(sv.i.rle,k=log(length(sv.i.rle$residuals))) # 1419
s2.sv.rle<-vcovHC(sv.i.rle)# "robust" SEs

## Crossvalidation exercise
cv.dat<-is.rs[rownames(sv.i.rle$model),]
rgrp<-y.hat.i<-y.hat.ni<-runif(nrow(cv.dat))
K<-10
for(j in 1:K){
  temp.out<-cv.dat[(rgrp>= (j-1)/K & rgrp<j/K),]
  temp<-cv.dat[setdiff(rownames(cv.dat), rownames(temp.out)),]
  fit.tmp0<-glm(sv.ni.rle$formula, data = temp, family = binomial(link = "probit"), 
    control = glm.control(epsilon = 1e-04, maxit = 100), model = T)
  fit.tmp1<-glm(sv.i.rle$formula, data = temp, family = binomial(link = "probit"), 
    control = glm.control(epsilon = 1e-04, maxit = 100), model = T)
  y.hat.ni[match(rownames(temp.out),rownames(cv.dat))]<-predict(fit.tmp0,newdata=temp.out,type="response")
  y.hat.i[match(rownames(temp.out),rownames(cv.dat))]<-predict(fit.tmp1,newdata=temp.out,type="response")
}

preds0<-prediction(y.hat.ni,cv.dat$svolik.dem)
preds1<-prediction(y.hat.i,cv.dat$svolik.dem)
perf0<-performance(preds0,measure="tpr", x.measure="fpr")
perf1<-performance(preds1,measure="tpr", x.measure="fpr")
plot(perf0, lwd=2)
plot(perf1,add=T, color="red", lty=3,lwd=2)
AuC0<-performance(preds0,measure="auc")#0.99
AuC1<-performance(preds1,measure="auc")#0.99

#getting results table
b.vars<-c("(Intercept)","log(labor.endow)","world.trade.gdp.interp", 
"log(labor.endow):world.trade.gdp.interp", "prop.dem.svolik", "nb.dem.prop.sv", 
"prev.dem.death.sv", "communist", "as.factor(xr.reg)goldStd","as.factor(xr.reg)interwar",
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.sv.m1 > 0)TRUE")

betas<-data.frame(coef(sv.i.rle)[b.vars])

row.names(betas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "Communist", "Gold std.", "Interwar",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.sv.rle[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(svolik.dem.m1)", "I(svolik.dem.m1 * log(labor.endow))",
"I(svolik.dem.m1 * world.trade.gdp.interp)", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)",
"I(svolik.dem.m1 * prop.dem.svolik)","I(svolik.dem.m1 * nb.dem.prop.sv)", 
"I(svolik.dem.m1 * prev.dem.death.sv)", "I(svolik.dem.m1 * gdppc.growth)")

gammas<-data.frame(coef(sv.i.rle)[g.vars])

row.names(gammas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.sv.rle[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.sv.rle["(Intercept)", "I(svolik.dem.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.sv.rle["log(labor.endow)", "I(svolik.dem.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.sv.rle["world.trade.gdp.interp", "I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.sv.rle["log(labor.endow):world.trade.gdp.interp", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["global prop. dem."])^2 + 2*s2.sv.rle["prop.dem.svolik","I(svolik.dem.m1 * prop.dem.svolik)"]
bgcov[6]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.sv.rle["nb.dem.prop.sv","I(svolik.dem.m1 * nb.dem.prop.sv)"]
bgcov[7]<-(betas.se["prior democratic failure"])^2 +2*s2.sv.rle["prev.dem.death.sv","I(svolik.dem.m1 * prev.dem.death.sv)"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model1betas.html",digits=3)
HTML(gammas,"model1gammas.html",digits=3)

## prep work for model interpretation

#scenarios
trade.span<-quantile(unique(is.rs$world.trade.gdp.interp),c(0.25,0.75),na.rm=T)
trade.span<-seq(trade.span[1],trade.span[2],length=150)
low.endow<-log(is.rs$labor.endow[is.rs$oid=="ARG1980"])
hi.endow<-log(is.rs$labor.endow[is.rs$oid=="CHN1980"])

low.endow.X<-cbind(1,0,low.endow,trade.span,0,0,
unique(is.rs$prop.dem.svolik[is.rs$year==1980]),0,
mean(is.rs$nb.dem.prop.sv[is.rs$year==1980]),0,0,0,0,0,0,1,0,0,
low.endow*trade.span,0)

hi.endow.X<-cbind(1,0,hi.endow,trade.span,0,0,
unique(is.rs$prop.dem.svolik[is.rs$year==1980]),0,
mean(is.rs$nb.dem.prop.sv[is.rs$year==1980]),0,0,0,0,0,0,1,0,0,
hi.endow*trade.span,0)

## effect of spatial terms
trade.med<-unique(is.rs$world.trade.gdp.interp[is.rs$year==1980])

loc.dem.lo<-cbind(1,0,mean(log(is.rs$labor.endow),na.rm=T),trade.med,0,0,
unique(is.rs$prop.dem.svolik[is.rs$year==1980]),0,
quantile(is.rs$nb.dem.prop.sv[is.rs$year==1980],.25),0,0,0,0,0,0,1,0,0,
mean(log(is.rs$labor.endow),na.rm=T)*trade.med,0)

loc.dem.hi<-cbind(1,0,mean(log(is.rs$labor.endow),na.rm=T),trade.med,0,0,
unique(is.rs$prop.dem.svolik[is.rs$year==1980]),0,
quantile(is.rs$nb.dem.prop.sv[is.rs$year==1980],.75),0,0,0,0,0,0,1,0,0,
mean(log(is.rs$labor.endow),na.rm=T)*trade.med,0)

pnorm(loc.dem.hi%*%coef(sv.i.rle))/pnorm(loc.dem.lo%*%coef(sv.i.rle))

no.loc.demtrans<-cbind(1,0,mean(log(is.rs$labor.endow),na.rm=T),trade.med,0,0,
unique(is.rs$prop.dem.svolik[is.rs$year==1980]),0,
quantile(is.rs$nb.dem.prop.sv[is.rs$year==1980],.5),0,0,0,0,0,0,1,0,0,
mean(log(is.rs$labor.endow),na.rm=T)*trade.med,0)

loc.demtrans<-cbind(1,0,mean(log(is.rs$labor.endow),na.rm=T),trade.med,0,0,
unique(is.rs$prop.dem.svolik[is.rs$year==1980]),0,
quantile(is.rs$nb.dem.prop.sv[is.rs$year==1980],.5),0,0,0,0,0,0,1,0,0,
mean(log(is.rs$labor.endow),na.rm=T)*trade.med,0)
loc.demtrans[17]<-1
pnorm(loc.demtrans%*%coef(sv.i.rle))/pnorm(no.loc.demtrans%*%coef(sv.i.rle))


#coefs

Betas<-mvrnorm(1000,coef(sv.i.rle),s2.sv.rle)
p.hat.low<-p.hat.hi<-NULL
for(i in 1:nrow(low.endow.X)){
	p.hat.low<-cbind(p.hat.low,pnorm(Betas%*%low.endow.X[i,]))
	p.hat.hi<-cbind(p.hat.hi,pnorm(Betas%*%hi.endow.X[i,]))
}

lo.mean<-apply(p.hat.low,2,mean)
hi.mean<-apply(p.hat.hi,2,mean)
lo.ci<-apply(p.hat.low,2,quantile,c(0.025,.975))
hi.ci<-apply(p.hat.hi,2,quantile,c(0.025,.975))

#Figure 4
#pdf("filename.pdf")
plot(trade.span,lo.mean,type="l",lwd=3,ylim=c(0,.08),
main="Predicted probabilities of democratic transitions",
ylab = "probability", xlab = "world trade, % of world GDP")
lines(trade.span,hi.mean,lty=2,lwd=3)
segments(trade.span,lo.ci[1,],trade.span,lo.ci[2,],lwd=1,col=gray(.6))
segments(trade.span,hi.ci[1,],trade.span,hi.ci[2,],lwd=1,col="blue",lty=3)
legend(x="topright",legend = c("labor scarce", "labor abundant" ), lty=c(1,2))
dev.off()

#model 2

sv.i.gdp<-glm(svolik.dem ~ I(svolik.dem.m1) +
              log(gdppc.use)*world.trade.gdp.interp +
              I(svolik.dem.m1*log(gdppc.use))*I(svolik.dem.m1*world.trade.gdp.interp) +
              prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) +
              nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) + prev.dem.death.sv +
              I(svolik.dem.m1*prev.dem.death.sv) + communist + as.factor(xr.reg) +
              I(nb.demtrans.sv.m1>0) + I(svolik.dem.m1*gdppc.growth), family=binomial(link="probit"),
              data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))
AIC(sv.i.gdp,k=log(length(sv.i.gdp$residuals)))
s2.sv.gdp<-vcovHC(sv.i.gdp)# "robust" SEs

# w/o RLE*trade interaction term for reference
sv.ni.gdp<-glm(svolik.dem ~ I(svolik.dem.m1) +
              log(gdppc.use)+world.trade.gdp.interp +
              I(svolik.dem.m1*log(gdppc.use))+I(svolik.dem.m1*world.trade.gdp.interp) +
              prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) +
              nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) + prev.dem.death.sv +
              I(svolik.dem.m1*prev.dem.death.sv) + communist + as.factor(xr.reg) +
              I(nb.demtrans.sv.m1>0) + I(svolik.dem.m1*gdppc.growth), family=binomial(link="probit"),
              data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))
AIC(sv.ni.gdp,k=log(length(sv.ni.gdp$residuals)))

#getting results table
b.vars<-c("(Intercept)","log(gdppc.use)","world.trade.gdp.interp", 
"log(gdppc.use):world.trade.gdp.interp", "prop.dem.svolik", "nb.dem.prop.sv", 
"prev.dem.death.sv", "communist", "as.factor(xr.reg)goldStd","as.factor(xr.reg)interwar",
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.sv.m1 > 0)TRUE")

betas<-data.frame(coef(sv.i.gdp)[b.vars])

row.names(betas)<-c("constant", "GDP", "world trade", "GDP*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "Communist", "Gold std.", "Interwar",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.sv.gdp[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(svolik.dem.m1)", "I(svolik.dem.m1 * log(gdppc.use))",
"I(svolik.dem.m1 * world.trade.gdp.interp)", "I(svolik.dem.m1 * log(gdppc.use)):I(svolik.dem.m1 * world.trade.gdp.interp)",
"I(svolik.dem.m1 * prop.dem.svolik)","I(svolik.dem.m1 * nb.dem.prop.sv)", 
"I(svolik.dem.m1 * prev.dem.death.sv)", "I(svolik.dem.m1 * gdppc.growth)")

gammas<-data.frame(coef(sv.i.gdp)[g.vars])

row.names(gammas)<-c("constant", "GDP", "world trade", "GDP*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.sv.gdp[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.sv.gdp["(Intercept)", "I(svolik.dem.m1)"]
bgcov[2]<-(betas.se["GDP"])^2 + 2*s2.sv.gdp["log(gdppc.use)", "I(svolik.dem.m1 * log(gdppc.use))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.sv.gdp["world.trade.gdp.interp", "I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["GDP*trade"])^2 + 2*s2.sv.gdp["log(gdppc.use):world.trade.gdp.interp", "I(svolik.dem.m1 * log(gdppc.use)):I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["global prop. dem."])^2 + 2*s2.sv.gdp["prop.dem.svolik","I(svolik.dem.m1 * prop.dem.svolik)"]
bgcov[6]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.sv.gdp["nb.dem.prop.sv","I(svolik.dem.m1 * nb.dem.prop.sv)"]
bgcov[7]<-(betas.se["prior democratic failure"])^2 +2*s2.sv.gdp["prev.dem.death.sv","I(svolik.dem.m1 * prev.dem.death.sv)"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model2betas.html",digits=3)
HTML(gammas,"model2gammas.html",digits=3)


#model 3
sv.i.ineq.rle1<-update(sv.i.rle,.~. + AW.90.10.interp + I(AW.90.10.interp^2) +
                       I(svolik.dem.m1*AW.90.10.interp) + I(svolik.dem.m1*(AW.90.10.interp^2)))
s2.sv.ineq.rle1<-vcovHC(sv.i.ineq.rle1)

# getting results table
b.vars<-c("(Intercept)","log(labor.endow)","world.trade.gdp.interp", 
"log(labor.endow):world.trade.gdp.interp", "prop.dem.svolik", "nb.dem.prop.sv", 
"prev.dem.death.sv", "AW.90.10.interp", "I(AW.90.10.interp^2)",
"communist", "as.factor(xr.reg)goldStd","as.factor(xr.reg)interwar",
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.sv.m1 > 0)TRUE")

betas<-data.frame(coef(sv.i.ineq.rle1)[b.vars])

row.names(betas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "ineq", "sq.ineq", "Communist", 
"Gold std.", "Interwar","post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.sv.ineq.rle1[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(svolik.dem.m1)", "I(svolik.dem.m1 * log(labor.endow))",
"I(svolik.dem.m1 * world.trade.gdp.interp)", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)",
"I(svolik.dem.m1 * prop.dem.svolik)","I(svolik.dem.m1 * nb.dem.prop.sv)", 
"I(svolik.dem.m1 * prev.dem.death.sv)", "I(svolik.dem.m1 * AW.90.10.interp)",                                           
"I(svolik.dem.m1 * (AW.90.10.interp^2))","I(svolik.dem.m1 * gdppc.growth)")

gammas<-data.frame(coef(sv.i.ineq.rle1)[g.vars])

row.names(gammas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "ineq", "sq.ineq", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.sv.ineq.rle1[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.sv.ineq.rle1["(Intercept)", "I(svolik.dem.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.sv.ineq.rle1["log(labor.endow)", "I(svolik.dem.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.sv.ineq.rle1["world.trade.gdp.interp", "I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.sv.ineq.rle1["log(labor.endow):world.trade.gdp.interp", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["global prop. dem."])^2 + 2*s2.sv.ineq.rle1["prop.dem.svolik","I(svolik.dem.m1 * prop.dem.svolik)"]
bgcov[6]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.sv.ineq.rle1["nb.dem.prop.sv","I(svolik.dem.m1 * nb.dem.prop.sv)"]
bgcov[7]<-(betas.se["prior democratic failure"])^2 +2*s2.sv.ineq.rle1["prev.dem.death.sv","I(svolik.dem.m1 * prev.dem.death.sv)"]
bgcov[8]<-(betas.se["ineq"])^2 +2*s2.sv.ineq.rle1["AW.90.10.interp","I(svolik.dem.m1 * AW.90.10.interp)"]
bgcov[9]<-(betas.se["sq.ineq"])^2 +2*s2.sv.ineq.rle1["I(AW.90.10.interp^2)","I(svolik.dem.m1 * (AW.90.10.interp^2))"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model3betas.html",digits=3)
HTML(gammas,"model3gammas.html",digits=3)



#model 4
sv.i.rle.ovb<-update(sv.i.rle, .~. + boix.relfract + boix.ethdiv + boix.prot + boix.cath
                     + boix.musl + I(svolik.dem.m1*boix.relfract) +
I(svolik.dem.m1*boix.ethdiv) + I(svolik.dem.m1*boix.prot) + I(svolik.dem.m1*boix.cath) + I(svolik.dem.m1*boix.musl))
s2.sv.rle.ovb<-vcovHC(sv.i.rle.ovb)

sv.i.ineq.rle2<-update(sv.i.rle.ovb,.~. + gini.sidd3+ I(gini.sidd3^2) + I(svolik.dem.m1*gini.sidd3) + I(svolik.dem.m1*(gini.sidd3^2)))
s2.sv.ineq.rle2<-vcovHC(sv.i.ineq.rle2)

# getting results table
b.vars<-c("(Intercept)","log(labor.endow)","world.trade.gdp.interp", 
"log(labor.endow):world.trade.gdp.interp", "prop.dem.svolik", "nb.dem.prop.sv", 
"prev.dem.death.sv", "gini.sidd3", "I(gini.sidd3^2)",
"communist", "as.factor(xr.reg)PostbretWood", "I(nb.demtrans.sv.m1 > 0)TRUE")

betas<-data.frame(coef(sv.i.ineq.rle2)[b.vars])

row.names(betas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "ineq", "sq.ineq", "Communist",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.sv.ineq.rle2[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(svolik.dem.m1)", "I(svolik.dem.m1 * log(labor.endow))",
"I(svolik.dem.m1 * world.trade.gdp.interp)", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)",
"I(svolik.dem.m1 * prop.dem.svolik)","I(svolik.dem.m1 * nb.dem.prop.sv)", 
"I(svolik.dem.m1 * prev.dem.death.sv)", "I(svolik.dem.m1 * gini.sidd3)",
"I(svolik.dem.m1 * (gini.sidd3^2))","I(svolik.dem.m1 * gdppc.growth)")

gammas<-data.frame(coef(sv.i.ineq.rle2)[g.vars])

row.names(gammas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "ineq", "sq.ineq", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.sv.ineq.rle2[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.sv.ineq.rle2["(Intercept)", "I(svolik.dem.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.sv.ineq.rle2["log(labor.endow)", "I(svolik.dem.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.sv.ineq.rle2["world.trade.gdp.interp", "I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.sv.ineq.rle2["log(labor.endow):world.trade.gdp.interp", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["global prop. dem."])^2 + 2*s2.sv.ineq.rle2["prop.dem.svolik","I(svolik.dem.m1 * prop.dem.svolik)"]
bgcov[6]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.sv.ineq.rle2["nb.dem.prop.sv","I(svolik.dem.m1 * nb.dem.prop.sv)"]
bgcov[7]<-(betas.se["prior democratic failure"])^2 +2*s2.sv.ineq.rle2["prev.dem.death.sv","I(svolik.dem.m1 * prev.dem.death.sv)"]
bgcov[8]<-(betas.se["ineq"])^2 +2*s2.sv.ineq.rle2["gini.sidd3","I(svolik.dem.m1 * gini.sidd3)"]
bgcov[9]<-(betas.se["sq.ineq"])^2 +2*s2.sv.ineq.rle2["I(gini.sidd3^2)","I(svolik.dem.m1 * (gini.sidd3^2))"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model4betas.html",digits=3)
HTML(gammas,"model4gammas.html",digits=3)


# Model 5
## DV = Polity

polity.i.rle<-glm(polity.democracy ~ I(polity.democracy.m1) + 
log(labor.endow)*world.trade.gdp.interp + I(polity.democracy.m1*log(labor.endow))*I(polity.democracy.m1*world.trade.gdp.interp) +
prop.dem.polity + I(polity.democracy.m1*prop.dem.polity) + 
nb.dem.prop.polity + I(polity.democracy.m1*nb.dem.prop.polity) + 
prev.dem.death.polity + I(polity.democracy.m1*prev.dem.death.polity) +
communist + as.factor(xr.reg) + I(nb.demtrans.polity.m1>0) + I(polity.democracy.m1*gdppc.growth), 
family=binomial(link="probit"), data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))
s2.polity.rle<-vcovHC(polity.i.rle)

# getting results table
b.vars<-c("(Intercept)","log(labor.endow)","world.trade.gdp.interp", 
"log(labor.endow):world.trade.gdp.interp", "prop.dem.polity", "nb.dem.prop.polity", 
"prev.dem.death.polity",
"communist", "as.factor(xr.reg)goldStd","as.factor(xr.reg)interwar",
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.polity.m1 > 0)TRUE")

betas<-data.frame(coef(polity.i.rle)[b.vars])

row.names(betas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "Communist", "Gold std.", "Interwar",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.polity.rle[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(polity.democracy.m1)", "I(polity.democracy.m1 * log(labor.endow))",
"I(polity.democracy.m1 * world.trade.gdp.interp)", "I(polity.democracy.m1 * log(labor.endow)):I(polity.democracy.m1 * world.trade.gdp.interp)",
"I(polity.democracy.m1 * prop.dem.polity)","I(polity.democracy.m1 * nb.dem.prop.polity)", 
"I(polity.democracy.m1 * prev.dem.death.polity)", "I(polity.democracy.m1 * gdppc.growth)")

gammas<-data.frame(coef(polity.i.rle)[g.vars])

row.names(gammas)<-c("constant", "RLE", "world trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.polity.rle[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.polity.rle["(Intercept)", "I(polity.democracy.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.polity.rle["log(labor.endow)", "I(polity.democracy.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.polity.rle["world.trade.gdp.interp", "I(polity.democracy.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.polity.rle["log(labor.endow):world.trade.gdp.interp", "I(polity.democracy.m1 * log(labor.endow)):I(polity.democracy.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["global prop. dem."])^2 + 2*s2.polity.rle["prop.dem.polity","I(polity.democracy.m1 * prop.dem.polity)"]
bgcov[6]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.polity.rle["nb.dem.prop.polity","I(polity.democracy.m1 * nb.dem.prop.polity)"]
bgcov[7]<-(betas.se["prior democratic failure"])^2 +2*s2.polity.rle["prev.dem.death.polity","I(polity.democracy.m1 * prev.dem.death.polity)"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model5betas.html",digits=3)
HTML(gammas,"model5gammas.html",digits=3)

# Model 6 w/ alternate trade.
## DV = Polity

polity.i.pwt<-glm(polity.democracy ~ I(polity.democracy.m1) + 
log(labor.endow)*EL.openk + I(polity.democracy.m1*log(labor.endow))*I(polity.democracy.m1*EL.openk) +
prop.dem.polity + I(polity.democracy.m1*prop.dem.polity) + 
nb.dem.prop.polity + I(polity.democracy.m1*nb.dem.prop.polity) + 
prev.dem.death.polity + I(polity.democracy.m1*prev.dem.death.polity) +
communist + as.factor(xr.reg) + I(nb.demtrans.polity.m1>0) + I(polity.democracy.m1*gdppc.growth), 
family=binomial(link="probit"), data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))
s2.polity.pwt<-vcovHC(polity.i.pwt)

# getting results table
b.vars<-c("(Intercept)","log(labor.endow)","EL.openk", 
"log(labor.endow):EL.openk", "prop.dem.polity", "nb.dem.prop.polity", 
"prev.dem.death.polity",
"communist", 
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.polity.m1 > 0)TRUE")

betas<-data.frame(coef(polity.i.pwt)[b.vars])

row.names(betas)<-c("constant", "RLE", "trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "Communist",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.polity.pwt[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(polity.democracy.m1)", "I(polity.democracy.m1 * log(labor.endow))",
"I(polity.democracy.m1 * EL.openk)", "I(polity.democracy.m1 * log(labor.endow)):I(polity.democracy.m1 * EL.openk)",
"I(polity.democracy.m1 * prop.dem.polity)","I(polity.democracy.m1 * nb.dem.prop.polity)", 
"I(polity.democracy.m1 * prev.dem.death.polity)", "I(polity.democracy.m1 * gdppc.growth)")

gammas<-data.frame(coef(polity.i.pwt)[g.vars])

row.names(gammas)<-c("constant", "RLE", "trade", "RLE*trade", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.polity.pwt[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.polity.pwt["(Intercept)", "I(polity.democracy.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.polity.pwt["log(labor.endow)", "I(polity.democracy.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["trade"])^2 + 2*s2.polity.pwt["EL.openk", "I(polity.democracy.m1 * EL.openk)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.polity.pwt["log(labor.endow):EL.openk", "I(polity.democracy.m1 * log(labor.endow)):I(polity.democracy.m1 * EL.openk)"]
bgcov[5]<-(betas.se["global prop. dem."])^2 + 2*s2.polity.pwt["prop.dem.polity","I(polity.democracy.m1 * prop.dem.polity)"]
bgcov[6]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.polity.pwt["nb.dem.prop.polity","I(polity.democracy.m1 * nb.dem.prop.polity)"]
bgcov[7]<-(betas.se["prior democratic failure"])^2 +2*s2.polity.pwt["prev.dem.death.polity","I(polity.democracy.m1 * prev.dem.death.polity)"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model6betas.html",digits=3)
HTML(gammas,"model6gammas.html",digits=3)


#### Checks:
#w/o hegemon. w/ RE. w/energy
sv.i.rle.noheg<-glm(svolik.dem ~ I(svolik.dem.m1) + 
log(labor.endow)*world.trade.gdp.interp + I(svolik.dem.m1*log(labor.endow))*I(svolik.dem.m1*world.trade.gdp.interp) +
prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) + 
nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) + 
prev.dem.death.sv + I(svolik.dem.m1*prev.dem.death.sv) +
communist + as.factor(xr.reg) + I(nb.demtrans.sv>0) + I(svolik.dem.m1*gdppc.growth), 
family=binomial(link="probit"), data= is.rs.noheg, model=T, control=glm.control(epsilon=0.0001, maxit=100))

sv.i.rle.re<-glmer(svolik.dem ~ I(svolik.dem.m1) + 
log(labor.endow)*world.trade.gdp.interp + I(svolik.dem.m1*log(labor.endow))*I(svolik.dem.m1*world.trade.gdp.interp) +
prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) + 
nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) + 
prev.dem.death.sv + I(svolik.dem.m1*prev.dem.death.sv) +
communist + as.factor(xr.reg) + I(nb.demtrans.sv>0) + I(svolik.dem.m1*gdppc.growth) + 
(1|WB.code), family=binomial(link="probit"), data= is.rs, model=T)


#using energy
sv.i.rle2<-glm(svolik.dem ~ I(svolik.dem.m1) + 
log(labor.endow)*world.trade.gdp.interp + I(svolik.dem.m1*log(labor.endow))*I(svolik.dem.m1*world.trade.gdp.interp) +
gw.energy2 + I(svolik.dem.m1*gw.energy2)+
prop.dem.svolik + I(svolik.dem.m1*prop.dem.svolik) + 
nb.dem.prop.sv + I(svolik.dem.m1*nb.dem.prop.sv) + 
prev.dem.death.sv + I(svolik.dem.m1*prev.dem.death.sv) +
communist + as.factor(xr.reg) + I(nb.demtrans.sv.m1>0) + I(svolik.dem.m1*gdppc.growth), 
family=binomial(link="probit"), data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))

s2.sv.rle2<-vcovHC(sv.i.rle2)# "robust" SEs
#getting results table
b.vars<-c("(Intercept)","log(labor.endow)","world.trade.gdp.interp", 
"log(labor.endow):world.trade.gdp.interp", "gw.energy2", "prop.dem.svolik", "nb.dem.prop.sv", 
"prev.dem.death.sv", "communist", "as.factor(xr.reg)goldStd","as.factor(xr.reg)interwar",
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.sv.m1 > 0)TRUE")

betas<-data.frame(coef(sv.i.rle2)[b.vars])

row.names(betas)<-c("constant", "RLE", "world trade", "RLE*trade", "energy", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "Communist", "Gold std.", "Interwar",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.sv.rle2[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(svolik.dem.m1)", "I(svolik.dem.m1 * log(labor.endow))",
"I(svolik.dem.m1 * world.trade.gdp.interp)", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)",
"I(svolik.dem.m1 * gw.energy2)", "I(svolik.dem.m1 * prop.dem.svolik)","I(svolik.dem.m1 * nb.dem.prop.sv)", 
"I(svolik.dem.m1 * prev.dem.death.sv)", "I(svolik.dem.m1 * gdppc.growth)")

gammas<-data.frame(coef(sv.i.rle2)[g.vars])

row.names(gammas)<-c("constant", "RLE", "world trade", "RLE*trade", "energy", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.sv.rle2[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.sv.rle2["(Intercept)", "I(svolik.dem.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.sv.rle2["log(labor.endow)", "I(svolik.dem.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.sv.rle2["world.trade.gdp.interp", "I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.sv.rle2["log(labor.endow):world.trade.gdp.interp", "I(svolik.dem.m1 * log(labor.endow)):I(svolik.dem.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["energy"])^2 + 2*s2.sv.rle2["gw.energy2","I(svolik.dem.m1 * gw.energy2)"]
bgcov[6]<-(betas.se["global prop. dem."])^2 + 2*s2.sv.rle2["prop.dem.svolik","I(svolik.dem.m1 * prop.dem.svolik)"]
bgcov[7]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.sv.rle2["nb.dem.prop.sv","I(svolik.dem.m1 * nb.dem.prop.sv)"]
bgcov[8]<-(betas.se["prior democratic failure"])^2 +2*s2.sv.rle2["prev.dem.death.sv","I(svolik.dem.m1 * prev.dem.death.sv)"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model3betas.html",digits=3)
HTML(gammas,"model3gammas.html",digits=3)

#using energy
# DV = polity

polity.i.rle2<-glm(polity.democracy ~ I(polity.democracy.m1) + 
log(labor.endow)*world.trade.gdp.interp + I(polity.democracy.m1*log(labor.endow))*I(polity.democracy.m1*world.trade.gdp.interp) +
gw.energy2 + I(polity.democracy.m1*gw.energy2) + prop.dem.polity + I(polity.democracy.m1*prop.dem.polity) + 
nb.dem.prop.polity + I(polity.democracy.m1*nb.dem.prop.polity) + 
prev.dem.death.polity + I(polity.democracy.m1*prev.dem.death.polity) +
communist + as.factor(xr.reg) + I(nb.demtrans.polity.m1>0) + I(polity.democracy.m1*gdppc.growth), 
family=binomial(link="probit"), data= is.rs, model=T, control=glm.control(epsilon=0.0001, maxit=100))
s2.polity.rle2<-vcovHC(polity.i.rle2)

# getting results table
b.vars<-c("(Intercept)","log(labor.endow)","world.trade.gdp.interp", 
"log(labor.endow):world.trade.gdp.interp", "gw.energy2", "prop.dem.polity", "nb.dem.prop.polity", 
"prev.dem.death.polity",
"communist", "as.factor(xr.reg)goldStd","as.factor(xr.reg)interwar",
"as.factor(xr.reg)PostbretWood", "I(nb.demtrans.polity.m1 > 0)TRUE")

betas<-data.frame(coef(polity.i.rle2)[b.vars])

row.names(betas)<-c("constant", "RLE", "world trade", "RLE*trade","energy", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "Communist", "Gold std.", "Interwar",
"post-Bretton Woods", "neigboring dem. transitions")

betas.se<-sqrt(diag(s2.polity.rle2[b.vars,b.vars]))
names(betas.se)<-row.names(betas)
betas[,2]<-betas.se
betas[,3]<-(1-pnorm(abs(betas[,1]/betas[,2])))*2

g.vars<-c("I(polity.democracy.m1)", "I(polity.democracy.m1 * log(labor.endow))",
"I(polity.democracy.m1 * world.trade.gdp.interp)", "I(polity.democracy.m1 * log(labor.endow)):I(polity.democracy.m1 * world.trade.gdp.interp)",
"I(polity.democracy.m1 * gw.energy2)", "I(polity.democracy.m1 * prop.dem.polity)","I(polity.democracy.m1 * nb.dem.prop.polity)", 
"I(polity.democracy.m1 * prev.dem.death.polity)", "I(polity.democracy.m1 * gdppc.growth)")

gammas<-data.frame(coef(polity.i.rle2)[g.vars])

row.names(gammas)<-c("constant", "RLE", "world trade", "RLE*trade", "energy", "global prop. dem.",
"prop. dem. neighbors", "prior democratic failure", "growth")

comb.vars<-row.names(gammas)[row.names(gammas)%in%row.names(betas)]

gammas[comb.vars,]<-gammas[comb.vars,] + betas[comb.vars,1]

gammas.se<-diag(s2.polity.rle2[g.vars,g.vars]) 
names(gammas.se)<-row.names(gammas)
bgcov<-NULL
bgcov[1]<-(betas.se["constant"])^2 + 2*s2.polity.rle2["(Intercept)", "I(polity.democracy.m1)"]
bgcov[2]<-(betas.se["RLE"])^2 + 2*s2.polity.rle2["log(labor.endow)", "I(polity.democracy.m1 * log(labor.endow))"]
bgcov[3]<-(betas.se["world trade"])^2 + 2*s2.polity.rle2["world.trade.gdp.interp", "I(polity.democracy.m1 * world.trade.gdp.interp)"]
bgcov[4]<-(betas.se["RLE*trade"])^2 + 2*s2.polity.rle2["log(labor.endow):world.trade.gdp.interp", "I(polity.democracy.m1 * log(labor.endow)):I(polity.democracy.m1 * world.trade.gdp.interp)"]
bgcov[5]<-(betas.se["energy"])^2 + 2*s2.polity.rle2["gw.energy2", "I(polity.democracy.m1 * gw.energy2)"]
bgcov[6]<-(betas.se["global prop. dem."])^2 + 2*s2.polity.rle2["prop.dem.polity","I(polity.democracy.m1 * prop.dem.polity)"]
bgcov[7]<-(betas.se["prop. dem. neighbors"])^2 + 2*s2.polity.rle2["nb.dem.prop.polity","I(polity.democracy.m1 * nb.dem.prop.polity)"]
bgcov[8]<-(betas.se["prior democratic failure"])^2 +2*s2.polity.rle2["prev.dem.death.polity","I(polity.democracy.m1 * prev.dem.death.polity)"]

gammas.se[comb.vars]<-gammas.se[comb.vars] + bgcov 
gammas[,2]<-sqrt(gammas.se)
gammas[,3]<-(1-pnorm(abs(gammas[,1]/gammas[,2])))*2

HTML(betas,"model7betas.html",digits=3)
HTML(gammas,"model7gammas.html",digits=3)

## END
